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ABSTRACT 

We present the final analysis of the data obtained 
by the Sydney University Giant Airshower Recorder 
(SUGAR). The data has been reanalysed to take into 
account the effects of afterpulsing in the 
photomultiplier tubes. Event data was used to 
produce a spectrum of "equivalent vertical muon 
number" and from this, a model dependent primary 
energy spectrum was obtained. These spectra show 
good evidence for the "Ankle": a flattening at 
1 0 1 9 eV. There is no sign of the cut-off 
which would be expected from the effects of the 
universal black body radiation. 

1 . Introduction 

The work was performed using an array which was operated at an 
atmospheric depth of 980 g cm on a site a few hundred km north 
of Sydney (latitude 30° 32 f south, longitude 149° 36’ east). 

The array and its results are described in some detail in (1) and (2). 

The geometry of the array was 54 points on a square grid where we 
operated autonomous "stations". At peak development, 47 of these 
points were occupied by operating stations. These were established on 
a set of nested square grids with spacings of 1600, 800 and 400 m. 

Each station had two liquid scintillator tanks buried 50 m apart in a 
north-south direction. The effective area of each scintillator was 6.0 
m 2 viewed by a single EMI 9623B photomultiplier tube. The 
threshold energy for detected muons was 0.75 sece GeV, where 6 is 
the zenith angle of the incident particles. The shape of each tank was 
designed so that a particle traversing any part of the scintillator at 
a given angle would produce the same light flux on the photomultiplier 
tube. 

The output (charge) pulse from the PM tube anode was deposited on a 
capacitor which then discharged with a decay time of 3. Ops. The 
potential difference across the capacitor was amplified and fed to a 
discr iminator set to trigger for a signal greater than or equal to 
that due to three coincident vertical muons passing through the 
scintillator. The width of the signal above the threshold was timed 
using a 10 MHz clock. Ideally, the width is proportional to the 
logarithm of the charge deposited by the PM tube. This arrangement, 
known as a logarithm height to time converter, is due to Suga and his 
coworkers. Its advantages are that it covers a wide dynamic range of 
signals and the number of clock pulses provides a convenient record of 
pulse size. The principal deficiency of this converter is that it 
makes the detector system prone to errors if the photomultiplier 
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afterpulses. This problem is further discussed below. 

The log height to time converter gave signals greater than the 
threshold at a rate of 30 per second; this was constantly monitored, 
as was the rate of pulses corresponding to more than 8 simultaneous 
vertical muons passing through the tank. 

Identical signal processing was used on both detector channels. If 
both discriminators at a station fired within 350 ns, a master trigger 
was generated and the following information was then recorded by a 
local tape recorder: 

(a) the time of the event as determined from the transmitted 
timing signal, 

(b) the widths of the pulses at the output of the discr iminators 
(recorded with 100 ns resolution). 

Such "local" events occured at a rate of about 12 per hour. At 
regular intervals simulated '’local' 1 events were generated by injecting 
simultaneous electrical pulses into the PM tube anode circuits of both 
detectors at the station. Four different sizes of injected pulses were 
used in sequence, covering the dynamic range of the system. 

2. Event Analysis 

The taped records of the local events were collected from all the 
stations, transferred to a computer and compared to find coincidences, 
within 80 \is t between three or more stations. Such coincidences 
are referred to as "array" events. Events were rejected from further 
analysis if the participating stations were collinear or if the signal 
times were not consistent with the passage of an air shower (the 
latter being classified "unphysical"). During the eleven years of 
operation 15327 events were recorded; of these 1238 were collinear, 

270 unphysical and 2 events were rejected on other grounds. Most of 
the events were registered in the minimum number of stations, namely 
three. 

The direction of the shower axis was obtained by the method of fast 
timing. If three stations participated in a shower, then a straight 
geometrical fit of a plane front was made to the timing data. If four 
or more stations participated in an array event, then a weighted least 
squares method was used to fit a plane shower front to the data. The 
number of muons in the shower and the location of the shower axis was 
found by a maximum likelihood method taking into account the records of 
all stations of the array (triggered or not). The likelihood 
expression consisted of: a term Q giving the probability of the 
observed responses of the triggered stations, a term P giving the 
probability of zero responses from non triggered stations and a term 
giving the a priori probability of the occurrence of a given shower 
size as proportional to raised to the power -3.0. 

The term Q was originally taken to be a Poissonian distribution with 
the mean set equal to the expected number of particles (above 100 



501 


particles this was replaced by a normal distribution with a fractional 
standard deviation of 10J). In the present analysis and that reported 
at the Bangalore Conference (3 and 4) we used a new formula for Q which 
was obtained after consideration of the possible effects of 
afterpulsing in the PM tubes of the array. Further details are given 
in the appendix. We feel quite confident that the new value of Q 
properly accounts for the effects of afterpulsing (5). 

The expected densities of muons were calculated using a structure of 
the Greisen type (6) modified to allow for the dependence of the outer 
slope parameter on zenith angle. 
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3. The Size Spectra 

2 til 

o.p < cos o g /.o y — T - 


The array was operated to produce a 
list of showers containing the 
following information: 

Time and date of event 
Zenith angle of shower axis 0 
Muon number N 
Position of aJcis X,Y 

In order to obtain shower spectra it 
is necessary to evaluate the array 
exposure; i.e. the effective area 
times the running time. We computed 
this using techniques substantially 
the same as described in Bell et al 
(1). A Monte Carlo method was used 
to simulate shower detection by our 
array; the probability of a station 
being operational was made to 
correspond to its average fractional 
on-time during an epoch. The total 
time of observation was divided into 
seven epochs during each of which 
the array configuration and station 
performance were substantially 
constant. The shower list and 
exposure information were combined 
to produce seven different spectra 
of N taken at different 
zenitfh angles, these are shown in 
figure 1. These spectra differ from 
one another because the threshold 
energy of the muons increases with 
the secant of the zenith angle and 
because showers of different zenith 
angle are at different stages of 
development, even though their 
primary energy may be the same. 

For the seven lowest N bins in each 
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Figure 1. The differential 
shower size (N ) spectra 
for various zeKith angle bins 
(6). The ordinate is the flux 
multiplied by N 3 . The 
lines are fitted by maximum 
likelihood to the lowest seven 
data points in each range of 0. 

of the zenith angle bands, 
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we used the maximum likelihood method to fit a differential power law 
spectrum of the form 

-V 0 

J[N u (e)] dN u - J e (N / N r ) dN y 
where N p * 3.16 x io e 

The five upper size bins were excluded from the fit since it appears 
that there may be a flattening of the primary energy spectrum for 
N > 10 8 . In any case, because of the small number of events 
iK these bins, the effect of including them in the fit is marginal. 

J and y were found to depend on cose in a linear 

manner allowing an analytic method (3) to be used to combine all the 
spectra into a single spectrum of the quantity N y : the equivalent 
vertical muon number, which is the size (N ) the shower 
would have had if it had entered the atmosphere vertically. 

The operation above is equivalent to using the well known ’’equal 
intensity cut method” (7) to obtain a shower size at any desired 
atmospheric depth. The spectrum of shown in figure 2 is 
for all showers with 0 < 60°. 



different zenith angles. Line (a) is the maximum likelihood fit to 
showers with log N < 7.75 and line (b) is a similar fit for log 
N > 7.75 Each poiPt represents the intensity in a (N ;0) 
bifn. Also shown is the result from the Akeno group y 
(1 < secQ < 1.1) and (c) is their fitted line. 

The line (a), with a differential slope of 3.35 ± 0.01 is fitted to 
all events with log N < 7.75. (If showers with 0 up to 
73° are included the ^lope steepens to 3-36.) 

It will be noticed that the equation of line (a) is effectively the 
same as presented by us in the 1 8th International Cosmic Ray 
Conference at Bangalore (3). Since then we have used an improved 
method of dealing with saturated detectors and this is the only 
difference between the two analyses. The effect of afterpulsing was 
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handled identically in the two cases. 

We also show on figure 2 the N spectrum from Akeno (8). We are 
able to explain the disagreement between our spectrum and the Akeno 
one by taking into account: 

(a) differing altitudes of the two arrays, 

(b) different threshold energies for entry of muons into the 
detectors , 

(c) different array geometries leading to muon densities being 
measured at different distances from the shower core, and 

(d) the use of different structure functions. 

In a similar way we are able to explain the disagreement between the 
our spectrum and those of the Nottingham and Yakutsk groups (9) and 
( 10 ). 

Note the flattening of the spectrum above log N = 8. A line 

fitted to events with log N >7.75 has a slope v of 

3.15 ± 0.14. One can compare the expected number of events which 

have log N > 8 according to line (a) (namely 58) with the 

observed ntimber (78). X i gives a probability of 1J that this 

should arise by chance. A stronger test, with the straight line fitted 

to all the data gives a x 2 probability of 3 %. 

. Energy Spectra 

The equivalent vertical muon size spectrum can be transformed into a 
primary energy spectrum using a conversion formula which is 
conventionally represented as a power law 

E = E r (N y /N r ) a 

N r is usually chosen to be in the middle of the range of shower 

sizes; in our calculations we used N =» 10 7 . 

r 

Two different models were used to obtain this conversion: the so- 

called Sydney model (11) and Hillas model E (12), In what follows we 
use the latter because of its general acceptance by workers in the 
field. 

The relevant relation is 
E = 1 .64 x 10 18 (N v /10 7 ) 1 •° 75 
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Figure 3 shows the spectrum produced. 
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Figure 3 Comparison of energy spectra from various groups and our 
spectrum converted by the Hillas E model (solid circles). Each 
spectrum is differential with the ordinate equal to the flux multiplied 
by E s . 

The squares are from the Leeds group at Haverah Park (13) 

The triangles are from the Utah group (14) 

The diamonds are from Volcano Ranch as quoted by (15) 

The inverted triangles are from the Yakutsk group (16) 

As with the N spectrum we have fitted two lines (which for clarity 
are not shown on figure 3), one below and one above 
log N = 7.75 with differential slopes of 
3.19 ¥ 0.01 and 2.99 ± 0.13 respectively. 

On figure 3 we also show the spectra from the Leeds group (13), 

Utah group (14), Volcano ranch as quoted by (15) and the Yakutsk group 
(16). The Leeds and the Volcano ranch spectra show the ankle 
feature and we confirm this. Another item of interest is the 
possibility of a spectrum cut-off at - 5 x 10 19 eV due to 
black body photons acting on particles from sources > 10 Mpc 
distant. There is no sign in our own spectrum of such a cut off. As 
described in section 3, the spectral slope flattens above 10 19 
eV and according to the Hillas E model we have eight showers above 
10 20 eV. 


5. Conclusions 

We have determined the spectrum of cosmic rays for energies 
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> 10 17 eV using data from a site at 31° south latitude. 

We used Hillas model E to find the energies of our events. 

The differential energy spectrum has a slope of 3 . 19 ± 0.01 below 
about 1 0 1 9 eV and the slope flattens to 2.99 ± 0. 1 3 above 
this energy. The spectrum extends beyond 10 20 eV with no sign 
of a cut off due to the Universal Black Body Radiation. 

Appendix - Afterpulsing 


An unfortunate feature of the logarithmic height to time converter is 
that pulses arriving late in a particular signal can keep it above 
threshold and hence cause the original signal height to be 
substantially overestimated. Ordinary random noise pulses from the PM 
tube are too infrequent to affect the results whereas pulses generated 
within a PM tube and occurring after a genuine air shower signal pulse 
can have serious effects. This phenomenon is known as afterpulsing. 

We looked at the possibility of excluding afterpulses by using 
electronic gating techniques but found that the inevitable gate 
pedestal created more problems than it solved. 


The seriousness of the effects of afterpulsing was not realised until 
the latter part of the experiment; even then, technology could not 
provide a practical solution. However we have developed a statistical 
method to account for the effects of afterpulsing. The technique 
incorporates the effects of afterpulsing and other signal enhancements 
into the probability distribution used in the shower fitting program. 
We describe below how this was done. 


During the last stages of the experiment we recorded for each station 
a sample of pulses, each of which had its height as well as it width 
measured. Samples were collected from all stations as part of the 
local event records and for a few stations, the comparison between 
pulse height and width was carried out for single tanks as well (17). 

From inspection of these data we chose a typical detector and removed 
its photomultiplier for further investigation of afterpulsing effects 
at larger pulse sizes. To this end we reactivated one of the original 
pilot arrray stations in Sydney and installed in it a PM tube with 
average afterpulsing. Data from this station was recorded for 8 days 
in a similar manner to the recordings made in the SUGAR array. To 
examine afterpulsing for pulses exceeding the height amplifier 
saturation level, the pulse height measuring system was connected to 
the output of the first stage of the preamplifiers rather than at the 
output of the final stage. We then ran the station for a further 
nine months in this "low gain" mode. The longer period was needed to 
accumulate an adequate data set at the lower rate. 

We used the results from the average tube as a model for all the 
others in the SUGAR array. An alternative approach would have been to 
establish a correlation between parameters derived from the data taken 
in the height/width runs and calibration data collected during the 
normal running of each station. If this had been successful; it would 
have allowed us to apply retrospective corrections to each station for 
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any time in the array 1 s operation. However we found no significant 
correlation allowing this to be done. Instead, the behaviour of the 
average tube was applied to all detectors of the array for the duration 
of the experiment. This was done by changing the probability 
distribution used in the shower analysis program. 

The data from the average tube was analysed statistically and a new 
formula was produced for the quantity Q used in the maximum likelihood 
expression (see main text). Q was previously Poissonian (with a mean 
set equal to the expected number of particles) for < 100 particles and 
was normally distributed with a fractional standard deviation of 
10? for > 100 particles. The new distribution for Q is of the gamma 
type. Its mean and standard deviation were fitted to the data from 
the average PM tube. The fit gave a constant fractional standard 
deviation of 30? and a mean which depended on the expected particle 
number n as 

(2.3 + 0.9n + 0.00 1 )n 2 )(1 .25 cose ) (0,i|0 “ 0,45 log n) 

where 0 is zenith angle of the particles (taken as equal to the 
zenith angle of the shower). 
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